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We investigate through molecular dynamics the transition from Knudsen to molecular diffusion 
transport towards 2d absorbing interfaces with irregular geometry. Our results indicate that the 
length of the active zone decreases continuously with density from the Knudsen to the molecular 
diffusion regime. In the limit where molecular diffusion dominates, we find that this length ap- 
proaches a constant value of the order of the system size, in agreement with theoretical predictions 
for Laplacian transport in irregular geometries. Finally, we show that all these features can be 
qualitatively described in terms of a simple random-walk model of the diffusion process. 
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The problem of Laplacian transport towards irregular 
surfaces represents an important subject of research in 
many fields of technological relevance including hetero- 
geneous catalysis, heat transfer and electrochemistry. In 
the particular case of catalysis, the role of the local sur- 
face morphology at the pore level on the global diffusion- 
reaction efficiency of the catalyst has been a subject of 
great interest in the past years [1—5]. Fractal concepts 
have been used to model complex surface roughness in 
the study of Knudsen diffusion in irregular reactive pores 
[6,7] including catalyst supports and adsorbents [8]. 

Another interesting point on the subject refers to the 
problem of the nonuniform accessibility of active sites 
located along an irregular reactive surface. If the sys- 
tem is diffusion-controlled, these so-called screening ef- 
fects may cause a dramatic reduction on the reactivity of 
the catalyst surface, as compared to the reactivity solely 
due to the intrinsic chemical reaction mechanism (i.e., 
the activity of the nominal surface). The extensive re- 
search developed on this field has been mainly devoted 
to the introduction, calculation and application of the 
concept of active zone in the Laplacian transport to and 
across irregular interfaces [9,10]. For example, through 
the coarse-graining method proposed in [9], it is possi- 
ble to determine the flux through an arbitrarily irregular 
surface from its geometry alone, avoiding the solution of 
the Laplace problem within a complex boundary domain. 
More recently, it has been shown that this technique pro- 
vides consistent predictions for the activity of catalyst 
surfaces [11]. 

All these studies dealing with active zone in Lapla- 
cian fields rely on the implicit assumption that molecu- 
lar diffusion is the governing mechanism of mass trans- 



port. Such an approximation, however, can only be lo- 
cally valid inside of the void space between the fins or ex- 
tended protrusions of an irregular surface if the mean free 
path of the diffusing molecules is sufficiently smaller than 
the width of these irregularities. For example, Knudsen 
diffusion may become the dominant mechanism of mass 
transport determining the reactivity of the system if the 
reagent is a diluted gas for which the collisions among 
molecules are less frequent than the collisions between 
the molecules and the catalytic surface [6,7]. The molec- 
ular mean free path therefore constitutes a lower cut- 
off for the validity of the molecular diffusion description. 
In a recent study [12], it has been shown analytically 
and through numerical simulations that surface rough- 
ness can have a significant effect on self-diffusion of gases 
in nanoporous media in the Knudsen regime. 

The aim of the present Letter is to study the transition 
in activity of an irregular absorbing interface when the 
mechanism of mass transport changes from Knudsen to 
molecular diffusion. Our approach is to use a nonequilib- 
rium molecular dynamics (NMD) technique in order to 
simulate a nonuniform and steady-state profile of reagent 
concentration between two active interfaces with an ar- 
bitrarily given roughness. This resembles more closely 
the experimental conditions frequently used in diffusion- 
reaction measurements. Here we adopt an NMD method 
that has been originally proposed for the study of self- 
diffusion in pure fluids [13]. The technique is entirely 
based on the standard molecular dynamics (MD) at equi- 
librium, but includes a special scheme to identify and 
exchange labeled and unlabeled particles during the sim- 
ulation. 

The MD part of the simulation consists in a two- 
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dimensional cell of size L x l x L y l containing N identical 
particles that interact through the Lennard- Jones poten- 
tial, $(Ary) = Ae[(a/Ar i3 ) 12 - (a/Ar i3 )% where Ar i3 
is the distance between particles i and j, e is the mini- 
mum energy, a is the zero of the potential, and we use 
I/a = 5.86. Periodic boundary conditions are applied 
in both the x and y directions. Distance, energy and 
time are measured in units of a, e and (ma 2 / 't) 1 / 2 , re- 
spectively, and the equations of motion are numerically 
integrated using the Verlct algorithm [14]. In all simula- 
tions we performed, the relative fluctuation around the 
average of the total energy of the system has always been 
smaller than 1CP 4 . 

After thermalization, two identical fractal interfaces 
of size Lyl and perimeter L p l are symmetrically placed 
into the system to simulate the roughness geometry of an 
absorbing material (see Fig. 1). At this point, the non- 
equilibrium dynamics is put forward through the follow- 
ing scheme: (1) half of the particles in the MD cell are 
randomly selected to carry a label, while the other half are 
left unlabeled] (2) every time a labeled (unlabeled) par- 
ticle crosses the interface at right (left) moving in the e" x 
(—e x ) direction it becomes unlabeled (labeled), and (3) 
when reinjected to the right (left) through the periodic 
boundaries in the x axis, an unlabeled (labeled) parti- 
cle becomes labeled (unlabeled). Although the particles 
are indistinguishable in terms of their mass and inter- 
action potential, this labeling technique naturally builds 
concentration gradients for both "species" that gradually 
develop to reach a desired non-uniform steady-state. In 
Fig. 1 we show the resulting stationary profiles along the 
x coordinate of the number fractions 6i = n;/(n; + n u ) 
and 9 U = 1 — 8i , where n; and n u are the number of la- 
beled and unlabeled particles, respectively, inside a ver- 
tical slice of fixed length in the system. From this point 
on during the simulation, we keep updating at each time 
step the number rii of particles being "absorbed" by the 
element i of the interface in order to compute its local 
mass flux ^ = rii/At, where At is the elapsed time after 
the steady-state has been established. Here we measure 
the efficiency of the interface in terms of the active length 
L a defined as [15] 

L„ = l/£>? (l<L a <L p ), (1) 

1=1 

where the sum is over the total number of interface ele- 
ments L p , and fa = qj is the normalized mass flux 
at element i. From the definition (1), L a = L p indicates 
a limiting state of equal partition of fluxes (fa = 1/L P , 
Vi) whereas L a of order L y should correspond to the case 
of an strongly "localized" flux distribution. 

Based on this NMD method, we performed simulations 
for 5 distinct initial configurations of the dynamical sys- 
tem, different values of the reduced temperature T, and 
reduced densities in the range 0.025 < p < 0.5, corre- 



sponding to systems with N — 1250 to 25000 particles. 
The evolution with time of the active length for any NMD 
realization reveals that, after a transient period, the sys- 
tem reaches a stationary state characterized by an aver- 
age value of L a that is representative of the flux distri- 
bution. As shown in Fig. 2, the active length decreases 
sharply with p for low density systems at T = 1.25, up 
to a point where it remains constant at L a w 27. The 
results from simulations performed at a higher temper- 
ature, T — 3.33, show that the behavior of the active 
length remains nearly the same, at least within the range 
of densities considered here. 

The decrease with density of the active length L a re- 
flects the transition from Knudsen to molecular diffusion 
in the distribution of activity at the interface. Because 
the mean free path of the particles for small p values 
is larger than the smaller length scale I of the irregu- 
lar interface, the activity is highly sensitive to geomet- 
rical constraints in the Knudsen regime. As shown in 
the inset of Fig. 2, the difference in behavior of the pair- 
correlation function for dilute and more dense fluids in- 
dicates that the simulated system undergoes noticeable 
structural modifications as p increases. 

At higher densities, the invariant behavior of L a is a 
consequence of molecular diffusion and can be explained 
in terms of an interesting theorem given by Makarov [16] 
to describe the properties of Laplacian fields on two- 
dimensional interfaces of arbitrary shape subjected to 
Dirichlet boundary conditions. The theorem states that 
the information dimension of the harmonic measure on a 
singly connected interface in d = 2 is exactly equal to 1. 
In terms of activity, this means that, regardless the shape 
of the interface, the total length L a of the region where 
most of the activity takes place should be of the order 
of the size L of the cell under a dilation transformation 
(see Ref. [15] for a detailed discussion of the active zone 
concept). Translating to our diffusion cell, where square 
Koch trees of third generation are the absorbing inter- 
faces, the theorem of Makarov predicts that the value of 
L a should be close to the size L v = 27, in good agreement 
with the NMD limit obtained for high densities. 

One can compute directly the length of the active zone 
for the two-dimensional continuum Laplacian problem 
which represents steady-state diffusion [17]. Numerical 
solutions of the Laplace equation, V 2 C = 0, for the con- 
centration field C inside the diffusion cell are obtained 
here through numerical discretization by means of finite- 
differences. For this, a constant unitary concentration is 
imposed at the source line (Co = 1) and Dirichlet bound- 
ary conditions (C = 0) are assigned to each elementary 
unit of the interface. Due to the symmetry with respect 
to the source line, only the concentration field in half of 
the domain needs to be calculated. The calculated value 
for the active length of the purely Laplacian cell is found 
to be L a = 22.9, compatible with the prediction of the 
Makarov theorem, L a ~ L y . 
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Now we show that a simple random-walk model of 
the diffusion-absorption process can provide a consistent 
qualitative description of the behavior observed in the 
NMD simulations. Adopting the same geometry, a parti- 
cle is released from a random position in the center line. 
The walker travels through the medium taking steps of 
random directions, but constant length A, till it crosses 
one of the wall elements of the fractal interface and gets 
absorbed. The flux at this element is then updated and 
the active length L a of the interface recalculated. For a 
fixed value of the step length A, the simulation goes on 
with particles being sequentially released and absorbed, 
till the active length reaches an average value that is ap- 
proximately constant. This value is usually obtained with 
less than 10 5 particles launched in the system. In Fig. 3 
we show the dependence on the parameter £ = (c/A) of 
the average L a computed for the third generation of the 
square Koch tree. For a two-dimensional gas, A can be 
interpreted as the mean free path, which is inversely pro- 
portional to the surface density of the system, A oc 1/ p. 
Similarly to the NMD simulations, two distinct regimes 
of activity can be clearly identified and directly related 
to the different governing mechanisms of mass transport, 
namely, Knudsen and molecular diffusion. At low values 
of £, the sharp decrease of L a reflects the strong influence 
on the mass transport process of the fractal geometry of 
the interface. The semi-log plot shown in the inset of 
Fig. 3 indicates that this decay in activity characteriz- 
ing the Knudsen regime is approximately logarithmic in 
shape. 

At sufficiently large values of £, the length L a reaches 
a plateau of minimum activity that is practically coin- 
cident with the value of the active length found for the 
Laplacian cell, L a — 22.9 (dashed line at the bottom 
in Fig. 3). It is important to show that these two ap- 
proaches to the problem provide consistent results for 
denser systems, even at the local scale of the interface 
geometry. Indeed, as displayed in Fig. 4, the normalized 
fluxes at each wall element generated from the continuum 
and random-walk methods are almost indistinguishable. 
Compared to the lower limit of the random-walk model, 
L a « 24, the higher value found for the active length 
with the NMD technique, L a « 27, can be explained in 
terms of the structural features of the simulated fluid. As 
shown in the inset of Fig. 2, the role of the attractive part 
of the interaction potential is to induce a smooth peak 
of short-range correlation that characterizes the presence 
of small particle clusters in the fluid. This geometrical 
aspect of the NMD system tends to increase the measure 
of the active length of the absorbing interface. 

In conclusion, we have shown through molecular dy- 
namics simulations that the active fraction of an irregu- 
lar absorbing interface should be sensitive to: (i) its ge- 
ometrical details; (ii) the governing mechanism of trans- 
port, and (hi) the structural aspects of the diffusing 
fluid. In particular, the active length for absorption 



of molecular diffusing fluids is found to be very close 
to that of a purely Laplacian system. These observa- 
tions may lead to new guidelines to the problem of dif- 
fusion and absorption on arbitrarily irregular interfaces. 
Furthermore, we have proposed a very simple random- 
walk model that incorporates the basic features of the 
diffusion-absorption process and is capable to describe, 
at least semi-quantitatively, the behavior of the active 
length for different diffusion regimes. This model pro- 
vides substantial insight on the effect of the diffusion 
mechanism on the interface activity and has the virtue of 
being computationally cheap. Finally, the approach in- 
troduced here is flexible enough to represent specific char- 
acteristics of irregular interfaces as well as other types of 
"absorption" mechanisms (e.g., finite-rate chemical reac- 
tions) limited by diffusion transport. As a consequence, 
the implications of this study are certainly relevant for 
the analysis and interpretation of diffusion and reaction 
processes occurring in a large variety of catalyst mate- 
rials. The modeling techniques we adopted should also 
be useful as design tools to choose a suitable catalytic 
interface for a given reaction-diffusion system. 
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FIG. 1. Schematic representation of the NMD diffusion 
cell. The absorbing interfaces are square Koch trees. Also 
shown in this figure is the dependence of the local number 
fraction 6 of labeled (full circles) and unlabeled (empty cir- 
cles) particles on the position along the x direction in the cell. 
In the case of the random-walk model, particles are released 
from random y positions at the dashed line in the center. 




FIG. 2. Dependence of the active length L a on the reduced 
density of the NMD cell for a fixed temperature, T = 1.25. 
The average values with error bars refer to simulations with 
5 different realizations of the NMD process. The horizontal 
dashed line at the top corresponds to the system size, L y = 27, 
while the one at the bottom indicates the value of the active 
length obtained from the simulation with the Laplacian cell, 
L a = 22.9. The inset shows the pair-correlation function cal- 
culated at the same temperature for a dilute (p — 0.025) and 
a more dense (p = 0.5) fluid. 
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FIG. 3. Dependence on the random-walk parameter £ of 
the active length L a of the fractal interface for the ran- 
dom-walk model. The horizontal dashed line at the top corre- 
sponds to the system size, L y = 27, while the one at the bot- 
tom gives the active length of the Laplacian cell, L a = 22.9. 
The semi-log plot in the inset shows that L a decays approx- 
imately logarithmically with £ in the regime of Knudsen dif- 
fusion. 
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FIG. 4. Distribution of the logarithm of the normalized 
fluxes crossing the wall elements along the absorbing frac- 
tal interface. The dashed line corresponds to a random-walk 
simulation performed with £ = 0.63, and the solid line is the 
distribution resulting from the numerical simulation with a 
Laplacian cell. 
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